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Abstract 

The  cumulant  theory  of  cyclosiationary  time-series 
is  applied  to  several  types  of  weak-signal  detection 
problems  that  arise  in  the  area  of  signal  interception, 
and  to  the  problem  of  estimating  the  relative  time- 
delay  of  a  heavily  corrupted  signal  that  is  received  at 
two  locations.  The  cumulant  theory  of  cyclostationar¬ 
ity  (CS)  is  the  theory  of  higher-order  temporal  and 
spectral  cumulants  and  moments  of  CS  time-series. 
Specifically,  the  theory  characterizes  the  additive  sine 
waves  present  in  the  output  of  nonlinear  transforma¬ 
tions  of  CS  time-series.  The  detection  and  time-delay 
estimation  problems  that  are  posed  are  difficult  to  solve 
because  the  signal  is  weak,  the  noise  and  interference  is 
nonstationary  and  non- Gaussian,  and  the  signal  does 
not  exhibit  second-order  CS. 

1  INTRODUCTION 

The  problems  considered  in  this  paper  can  be  for¬ 
mulated  in  terms  of  the  following  two-sensor  signal 
model 

M 

x{t)  =  s(t)-fm,(f)-b53i*(f)  (1) 

f  k=\ 

>  M 

y{t)  =  Aos(t-|-do)  +  ini,(t)  +  ^Afcifc(t-|-(it), 

t=i 

where  s(t)  is  the  signal  of  interest  (SOI),  m,(f)  and 
my(f)  are  independent  white  Gaussian  noises,  the 
are  signals  not  of  interest  (SNOIs),  {>!*}  are 
attenuation  factors,  and  {</*}  are  the  relative  delays 
between  the  signal  components  in  x{t)  and  y{t).  For 
example,  suppose  that  M  =  2  and  the  power  levels 
of  m,  and  are  time-varying.  Further,  assume  that 
the  SOI  is  weak  and  fourth-order  cyclostationary  (CS) 

*This  work  was  jointly  supported  by  the  National  Science 
Foundation  under  grant  MIP-91-12800  and  the  Army  Research 
Office  under  contract  DAAL03-91-C-0018. 


-  but  not  second-order  CS  [8]  -  with  period  To,  and 
the  two  interferers  are  second-order  CS  with  periods 
Ti  and  Tj  such  that  Tj  +  r2  =  Tq.  The^rst  prob¬ 
lem  is  to  detect  the  presence  of  the  SOI  given  a  finite 
segment  of  x(t).  The  next  problem  is  to  determine 
the  parameter  do  given  finite  segments  of  both  x{t) 
and  y{t).  These  problems  are  difiScult  to  solve  using 
the  stationary  models  of  the  various  signals  (which 
leads  to  radiometry  [4]  and  generalized  cross  correla¬ 
tion  methods  [6])  because  for  detection  the  nonstation¬ 
ary  noise  and  interference  complicates  the  threshold 
setting,  and  for  time-delay  estimation  the  interference 
corrupts  the  relevant  phase  information  in  the  cross 
spectrum.  The  theory  of  higher-order  statistics  is  not 
helpful  because  it  also  is  not  signal  selective:  all  of 
the  signals  contribute  to  higher-order  cumulants  for 
a  stationary  signal  model  [8J.  The  theory  of  second- 
order  CS  (SOCS)  is  also  not  helpful  because  the  SOI 
has  no  second-order  cyclic  features  [3].  However,  these 
two  tasks  (and  other  similar  tasks)  can  be  handled  by 
using  the  higher-order  CS  (HOCS)  of  the  SOI.  Before 
showing  how  this  can  be  done,  we  first  review  the  defi¬ 
nitions  of  higher-order  moments  and  cumulants  for  CS 
time-series. 

2  THE  PARAMETERS  OF  HOCS 

For  the  time-series  x{t)  for  -oo  <  t  <  oo,  we  define 
the  nth-order  lag-product  time-series  by 

Lz{t,  t)„  =  J][  a;(t  -b  Tj),  T  =  [n  •  •  •  r„]^  (2) 

>=i 

where  f  denotes  transposition.  The  nth-order  cyclic 
temporal  moment  function  (CTMF)  is  defined  by  the 
limiting  time  average 


KMn  =  ^  /  Lz{t,T)n 

T—oo  1 

=  <L,(f,T)„e-«-‘)„, 


dt 


and  is  simply  the  Fourier  coeflScient  associated  with 
the  component  in  the  time-series  r)„.  It 

can  be  seen  that  the  CTMF  is  a  Fourier  coefficient 
of  a  moment  function  because  the  nth-order  fraction- 
of-time  probabilistic  moment  (the  temporal  moment 
function)  associated  with  Lx{t,T)n  can  be  expressed 
as  [1] 

=  (4) 

O' 

where  the  sum  is  over  all  real  numbers  a,  called  nth- 
order  impure  cycle  frequencies,  for  which  R^(t)„  ^  0. 
The  functions  (3)  and  (4)  exist  and  are  well-behaved 
for  appropriate  models  of  many  time-series  including 
amplitude  modulated,  pulse-amplitude-modulated, 
and  phase-shift-keyed  and  frequency-shift-keyed  sig¬ 
nals  [8]. 

The  nth-order  temporal  cumulant  function  (TCF) 
for  the  set  of  time-series  translates  {x(t  -f-  is 

defined  by  (cf.  [7]  and  [8]) 

p 

Cx{t,T)n=  ^  ^  h(p)'j^Rx(t,T,,j)nj  .  (5) 

L  •’■=1 

The  sum  in  (5)  is  over  the  set  P  of  all  distinct  parti¬ 
tions  of  the  index  set  {l,2,---,n},  where  each  par¬ 
tition  has  p  elements,  1  <  p  <  n,  and 

^(p)  —  1)”  The  vector  t„.  consists  of  the 

Uj  lags  with  indices  in  the  set  uj.  The  cyclic  temporal 
cumulant  function  (CTCF)  is  the  Fourier  coefficient 
of  the  TCF; 

=  (6) 

where  /?  is  an  nth-order  pure  cycle  frequency  for  x(t) 
[2,  8].  Combining  (4)-(6)  reveals  that  the  CTCF  is 
given  by  the  following  explicit  function  of  lower-order 
CTMFs: 


for  many  time-series  of  interest  and,  therefore,  it  is 
strictly  Fourier  transformable  [8].  The  RD-CTCF 
is  simply  the  CTCF  associated  with  the  n  variables 
{x(t  -I-  with  r„  =  0.  We  use  the  notation 

Cf  («)„  =  Cf  (r)„  r  =  [ti0],  (7) 

where  u  is  an  (n-l)-dimensional  vector.  The  (n-1)- 
dimensional  Fourier  transform  of  (7)  is  denoted  by 

/oo 

Cf(«)„e-*2'«'/  du,  (8) 

•OO 

where /'  =  [/i  ••• 

Consider  the  n  complex-demodulate  time-series 
XT{tifj)  for  j  =  l,'”,n,  associated  with  narrow 
bandpass  filtered  versions  of  x(f),  where 

ft+Tl2 

XT{t,f)=  xiv)e-'^^f'’ dv.  (9) 

Jt-TI2 

The  limit  as  T  — >  (X>  of  the  limiting  time-average  of 
the  product  of  these  spectral  components  is  called  the 
spectral  moment  function 

Sxif)„^l^^(jlXTit,fj)'j  ,  (10) 

and  it  can  be  shown  that  Dirac  deltas  (denoted  by 
6(*))  can  be  factored  out  as  follows: 

5»(/)„  =  J^5“(/%5(/^i-a).  (11) 

or 

However,  the  factor  5“(/0n  typically  contains  addi¬ 
tional  Dirac  deltas  [8].  The  spectral  cumulant  function 
is  given  by 


Cf(r)„=  *(P)  E 

^=^1=1  an=0j=i 

where  i  =  [1---1]^,  and  a  =  [al•••ap]^  An  nth- 
order  pure  cycle  frequency,  as  first  defined  in  [2],  is 
the  frequency  of  a  finite-strength  additive  sine-wave 
component  in  i*(t,  t)„  that  is  free  of  all  contributions 
from  products  of  sine-wave  components  from  lower- 
order  lag  products  obtained  by  factoring  £*(f,r)„. 

The  CTMF  and  the  CTCF  are  not  in  general  ab¬ 
solutely  integrable  due  to  the  presence  of  sinusoidal 
components  in  t.  These  components  formally  result 
in  Dirac  deltas  in  the  n-dimensional  Fourier  trans¬ 
form  of  the  CTCF.  However,  a  reduced-dimension 
(RD)  version  of  the  CTCF  is  absolutely  integrable 


Px(/)„=  E  %)n5x(/.,)n,  ,  (12) 

where  is  the  vector  of  frequencies  with  subscripts 
in  the  set  vj,  and  it  follows  from  (8)  and  (11)  that 
Dirac  deltas  can  be  factored  out: 

Px(/)„  =  E  (/')n«(/^l  -  /9),  (13) 

where  Px{f')n  is  called  the  cyclic  polyspectrum  (CP) 
and  does  not  contain  Dirac  deltas.  As  first  shown  in 
[2],  the  CP  is  the  (n-l)-dimensional  Fourier  transform 
of  the  RD-CTCF  (8).  This  is  a  generalization  of  the 
Wiener  relation  between  the  power  spectrum  and  au¬ 
tocorrelation  from  second-order  stationary  time-series 


to  nth-order  CS  time-series  (cf.  [1]).  For  strictly  sta¬ 
tionary  time-series,  the  only  cycle  frequency  (pure  or 
impure)  is  a  =  =  0.  For  a  strictly  stationary 

stochastic  process,  the  same  is  true  for  probabilistic 
versions  of  the  temporal  moments  and  cumulants  de¬ 
fined  in  this  section,  which  form  the  basis  for  the  the¬ 
ory  of  higher-order  statistics.  Such  parameters  are  not 
signal  selective,  as  noted  in  the  next  section,  whereas 
the  cyclic  parameters  (/3  7^  0,  a  ^  0),  which  form  the 
basis  for  the  theory  of  HOCS,  can  be  signal  selective. 

3  SIGNAL  SELECTIVITY 


An  important  advantage  of  exploiting  CS  in  sig¬ 
nal  processing  tasks  is  that  the  cyclic  parameters  are 
signal  selective  in  that  the  parameters  associated  with 
the  SOI  can  be  estimated  from  data  that  also  contains 
noise  and  SNOIs;  and  as  the  amount  of  data  becomes 
infinite,  the  effects  on  the  estimate  of  the  noise  and 
interference  vanish.  The  nature  of  the  signal  selec¬ 
tivity  properties  of  higher-order  cyclic  moments  and 
cumulants  is  examined  next. 

Because  the  signals  and  noises  in  (1)  are  assumed 
to  be  mutually  independent,  the  TCF  for  x{t)  is  given 
by  the  simple  formula  [8] 

M 

C'*(t,  T)n  —  Cg{t,  T)n  -b  Cm,  {t,T)n  +  ^  \  t)„, 

k=l 

which  implies  that  if  /3  is  an  nth-order  pure  cycle  fre¬ 
quency  for  s(t),  and  is  not  for  any  of  the  other  signals 
(it  is  unique  to  s(f)),  then  C^(t)„  =  Simi¬ 

larly,  if  Pk  is  a  unique  nth-order  pure  cycle  frequency 
for  ik{t),  then  =  Cf{T)„. 

For  moments,  the  signal  selectivity  property  de¬ 
pends  on  the  cycle  frequencies  for  all  of  the  signals 
present  for  all  orders  m  <  n.  This  can  be  seen  by 
expressing  the  CTMF  in  terms  of  CTCFs  [8]: 


E 


p 


E 


■  (14) 


If  there  is  a  vector  /3  such  that  at  least  one  of  its  el¬ 
ements  is  a  cycle  frequency  that  is  associated  with  a 
SNOI,  then  the  CTMF  for  x{t)  will  not  be  equal  to 
the  CTMF  for  s(t).  Nevertheless,  it  can  happen  that 
the  contributions  to  the  CTMF  from  SNOIs  do  not 
affect  the  phase  of  the  CTMF  but  only  its  magnitude. 
For  moments,  then,  there  are  two  kinds  of  signal  selec¬ 
tivity,  depending  on  what  information  in  the  moment 
is  considered  useful  (magnitude  or  phase).  Because  of 
this  potentially  troublesome  complication,  the  signal- 
selectivity  properties  of  cumulants  are  deemed  more 


useful  than  those  for  moments.  Thus,  the  following 
sections  describe  algorithms  in  terms  of  cumulants, 
but  in  most  cases,  an  alternative  algorithm  can  be 
created  by  simply  replacing  cumulants  with  moments. 
The  usefulness  of  the  resulting  alternatives  depends 
on  the  cycle  frequencies  associated  with  the  SNOIs, 
which  must  be  evaluated  on  a  case-by-case  basis. 

4  WEAK-SIGNAL  DETECTION 

In  this  section,  the  problem  of  detecting  the  pres¬ 
ence  of  the  signal  s{t)  in  a  received  data  set  x{t)  as  in 
(1)  is  considered.  There  are  several  versions  of  this  de¬ 
tection  problem  that  are  of  interest.  The  first  is  called 
the  general  search  problem,  in  which  a  data  set  is  ana¬ 
lyzed  to  determine  if  there  are  any  CS  signals  present. 
No  information  about  the  received  data  is  assumed  to 
be  known  in  the  genercd  search  problem.  In  the  second 
problem,  called  the  knoum- cycle- frequency  problem,  a 
specific  pure  cycle  frequency/order  pair  {p,no)  is  of 
interest,  and  it  is  desired  to  determine  if  there  is  a  sig¬ 
nal  present  in  the  data  corresponding  to  this  pair.  In 
the  third  problem,  called  the  known-modulation  prob¬ 
lem,  the  modulation  format  of  the  signal  of  interest  is 
known,  and  hence  the  cyclic  cumulants  of  the  signal 
are  known  (in  principle);  it  is  desired  to  determine  the 
presence  or  absence  of  this  signal. 

The  General  Search  Problem 
In  this  problem  there  is  a  maximum  order  N  of  non¬ 
linearity  that  is  to  be  used  for  processing.  The  goal 
of  the  processing  is  to  produce  a  list  of  pure  cycle  fre¬ 
quencies  {/?«}  for  each  order  n  <  N.  This  list  {Pn} 
for  each  n  characterizes  the  detectable  CS  of  order  n 
(and  only  n)  that  is  associated  with  x(t).  Thus,  these 
lists  are  not  contaminated  by  entries  that  are  due  to 
lower-order  sine  wave  interactions.  To  accomplish  this 
task,  the  TCF  is  estimated  for  x(t)  for  each  order  n. 
From  the  estimate  of  the  TCF  of  order  n,  the  cycle 
frequencies  {Pn},  which  are  needed  for  the  estimate  of 
the  TCF  for  order  n  -f  1,  can  be  found.  More  explic¬ 
itly,  the  general  search  problem  can  be  tackled  using 
the  following  algorithm: 

0.  Let  71  =  1 

1-  ^x(^>  ‘^)n  ~  ‘*’)n  ~  ^  ]  Cx{t,  T )„j 

P 

2.  r(/)  =  FFT,  {c;(t,T)„} 

3.  Threshold  detect  the  bins  of  Y  to  find  {/?„} 

4.  Cf"(T)„  =  (ci(t,r)„e-*2'/’»‘)^ 

5.  C'*(t,r)„  =  ^Cf»(r)„e'2-^»‘ 

6.  71  — >  n  -I- 1;  if  71  <  iV  then  go  to  1. 


In  step  4,  the  interval  over  which  the  average  {•)q  is 
performed  is  determined  by  the  amount  of  data  avail¬ 
able.  If  any  of  the  detected  cycle  frequencies  are  of 
particular  interest,  a  cyclic  polyspectral  analysis  can 
be  performed  from  which  the  modulation  type  can 
possibly  be  determined  [8]. 

The  Known- Cycle- Frequency  Problem 
In  this  problem,  one  or  more  of  the  signal’s  modulation 
frequencies,  such  as  a  symbol  rate  or  carrier  frequency, 
is  assumed  to  be  known,  but  the  exact  functional  form 
of  the  CTCF  is  unknown.  The  environment  is  still 
assumed  to  be  unknown  and,  therefore,  the  general 
search  algorithm  is  still  of  interest.  However,  it  can 
be  improved  for  the  known-cycle-frequency  problem 
by  combining  it  with  a  least-squares  estimation  tech¬ 
nique.  Let  {P,  no)  be  the  cycle  frequency /order  pair  of 
interest.  Use  the  general  search  algorithm  up  to  order 
no  —  1-  Form  C''(f,  t)„„,  and  use  a  least-squares  esti¬ 
mator  to  detect  the  presence  of  the  signal  of  interest 
using  the  statistic 

where 

and  where  w  is  the  least-squares  weight  vector 

w  =  arg  imn  ^  .  (15) 

The  solution  to  (15)  is  th  =  (r)„„,  where 

il=(cl(f,T)„„cl(f,r)")^, 

in  which  H  denotes  conjugate  transpose.  Thus,  the 
detection  statistic  is 

r  =  cf(r)"Jl-^cf(r)„„, 

which  is  obtained  by  forming  the  particular  linear 
combination  of  data  sets  C  (f,  ri)„„ ,  •  •  • ,  C  (t, TK)„g 
that  optimally  combines  the  regenerated  sine  waves 
with  frequency  /?  present  in  each  set,  and  then  cor¬ 
relates  this  composite  regenerated  sine  wave  with  the 
stored  sine  wave 

The  Known-Modulation  Problem 
In  this  problem,  it  is  desired  to  determine  if  a  signal 
with  known  modulation  type  is  present.  In  particular, 
the  CTCF  of  s(t)  for  n  =  uq  and  pure  cycle  frequency 
p  is  known.  The  general  search  algorithm  can  be  used 


to  remove  all  lower-order  sine  waves  up  to  order  no  —  1. 
Then,  from  C'^{t,r)„^  the  CTCF  estimate 
for  cycle  frequency  P  can  be  determined  by  computing 
the  Fourier  coeflScient  as  in  (6).  The  proposed  detec¬ 
tion  statistic  is 

y=  r 

J  — OO 

The  primary  justification  for  this  particular  statistic 
is  that  when  no  signal  is  present  with  nth-order  pure 
cycle  frequency  P,  then  0,  which  implies 

that  y  — >  0;  when  the  signal  of  interest  is  present, 
then 

/OO 

|Cf(u)„|"  du.  (16) 

•OO 

Thus,  j/  is  an  asymptotically  noise-free  statistic  on 
both  the  signal-present  and  signal-absent  hypotheses. 
Furthermore,  the  integral  (16)  is  finite  [8].  Hence,  this 
statistic  is  the  natural  generalization  of  the  single  cy¬ 
cle  detector  that  exploits  SOCS  [4]. 

The  detection  statistic  Y  can  be  generalized  to  in¬ 
clude  only  a  portion  of  «-space,  denoted  by  G  C  R”° , 

y  =  /  c^Au)„,c^{u)*„ju. 
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Choices  for  G  might  include  those  values  of  u  for  which 
the  RD-CTCF  («)„„  is  particularly  large,  or  for 
which  the  coeflBcient  of  variation  (variance  divided  by 
squared  mean)  of  the  estimator  (7^j.(u)„„  of  the  RD- 
CTCF  is  particularly  small  [8]. 

5  TIME-DELAY  ESTIMATION 

Conventional  approaches  to  the  problem  of  esti¬ 
mating  the  time-delay  (or  time  difference  of  arrival 
(TDOA))  between  signal  components  in  data  from 
two  sensors  can  be  collectively  referred  to  as  gener¬ 
alized  cross  correlation  (GCC)  methods  [6].  In  the 
GCC  methods,  filtered  versions  of  the  sensor  outputs 
x{t)  and  y{t)  are  cross  correlated,  and  the  estimate 
of  do  is  taken  to  be  the  location  of  the  peak  in  the 
cross-correlation  estimate.  These  methods  suffer  when 
interferers  are  present  (M  >  1  in  (1)),  because  each 
interferer  contributes  a  peak  of  its  own  to  the  cross  cor¬ 
relation  function.  This  causes  two  problems.  The  first 
is  a  resolution  problem  which,  to  be  solved,  requires 
that  the  differences  in  the  TDOAs  for  each  of  the  sig¬ 
nals  be  greater  than  the  widths  of  the  cross  correla¬ 
tion  functions  so  that  the  peaks  can  be  resolved.  The 
second  problem  is  that  it  is  difficult  to  correctly  as¬ 
sociate  each  peak  with  its  corresponding  signal.  Both 
of  these  problems  arise  because  the  GCC  methods  are 
not  signal  selective;  they  produce  TDOA  peaks  for 


all  the  signals  in  the  received  data  unless  they  are 
spectrally  disjoint  and  can,  therefore,  be  separated 
by  filtering.  Signal-selective  methods  that  exploit  the 
SOCS  of  the  desired  signal,  which  is  assumed  to  be 
unique  to  that  signal,  are  studied  in  [5].  These  meth¬ 
ods  have  been  shown  to  outperform  the  GCC  methods, 
and  have  been  shown  to  produce  unbiased  TDOA  esti¬ 
mates  with  variance  that  is  smaller  than  the  Cramer- 
Rao  lower  bound  on  the  variance  of  TDOA  estimators 
that  are  based  on  the  assumption  that  the  signal  and 
its  environment  are  stationary.  However,  these  meth¬ 
ods  fail  when  there  is  no  SOCS  to  exploit.  In  this  case, 
the  theory  of  HOCS  can  be  used  to  develop  signal- 
selective  TDOA  estimators.  Following  the  approach 
in  [5]  for  SOCS,  the  methodology  considered  here  for 
HOCS  is  based  on  least-squares  estimation.  The  fol¬ 
lowing  two  examples  illustrate  the  methodology  [8]. 

Define  a  cross  cumulant  between  n-1  time- 
translates  of  x{t)  and  one  translate  of  y(t)  as  follows 

Cxyit,  t)„  =  Cumulant  {y(t  -f  r„),  x{t  +  rj)}"=i 


which  suggests  the  following  least-squares  approach: 

do  =  arg  min  f  \C^^{u)„  -  du. 
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The  estimator  for  do  is  given  explicitly  by 

do  =  (17) 

which  is  a  higher-order  generalization  of  the  second- 
order  Cyclic  Phase  Difference  algorithm  for  TDOA  es¬ 
timation  without  cross-sensor  measurements  [5]. 

6  CONCLUSIONS 

The  higher-order  cumulants  and  moments  of  cyclo¬ 
stationary  time-series  can  be  used  to  perform  diffi¬ 
cult  detection  and  estimation  tasks.  The  basic  signal- 
selectivity  properties  of  cumulants  yield  estimators 
that  are  tolerant  to  noise  and  cyclostationary  inter¬ 
ference.  This  admittedly  terse  presentation  will  be 
expanded  on  in  a  forthcoming  journal  paper  (see  [8]). 


The  Fourier  coefficient  of  this  cross  cumulant  for  the 
cycle  frequency  /3  for  the  signal  model  (1)  (assuming 
that  the  noise  and  interference  do  not  exhibit  nth- 
order  CS  with  pure  cycle  frequency  /?)  is  given  by 

Ci(r)„  ^ 

=  AoC^{T  +  6ndo)„, 

where  S„  is  the  unit  vector  along  the  nth  coordinate. 
It  is  easy  to  show  that  the  following  relations  involving 
RD-CTCFs  hold: 

C^yiu)„  =  AoCf («-ido)ne‘=''^‘'‘’;  Cf («)„  =  Cf(«)„ 

This  suggests  a  least-squares  fit  of  a  measurement  of 
Cxy  to  a  measurement  of  over  a  region  G  of  u- 
space  of  interest: 


mm 

A 


“/ 

A  Jg 


du, 


which  leads  to  the  following  estimator  of  the  delay  do¬ 


do  =  arg  max  3?  [  Cf  (tt)„C^  (« -f  id)* d«. 

a  Jg 


This  estimator  is  a  higher-order  generalization  of  the 
SPECtral  Coherence  Alignment  algorithm  for  TDOA 
estimation  [5],  which  exploits  SOCS,  and  has  been 
shown  to  possess  several  optimality  properties. 

As  an  alternative,  cross-sensor  measurements  can 
be  avoided  entirely  by  noting  that 

C^iu)n  =  («)„e'2"^‘'»;  cf(«)„  =  cf  («)„, 
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